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The relaxation process of one dimensional surface modulations is re- 
examined. Surface evolution is described in terms of a standard step flow 
model. Numerical evidence that the surface slope, D(x,t), obeys the scaling 
ansatz D(x, t) = a(t)F(x) is provided. We use the scaling ansatz to transform 
the discrete step model into a continuum model for surface dynamics. The 
model consists of differential equations for the functions a(t) and F(x). The 
solutions of these equations agree with simulation results of the discrete step 
model. We identify two types of possible scaling solutions. Solutions of the 
first type have facets at the extremum points, while in solutions of the second 
type the facets are replaced by cusps. Interactions between steps of opposite 
signs determine whether a system is of the first or second type. Finally, we 
relate our model to an actual experiment and find good agreement between a 
measured AFM snapshot and a solution of our continuum model. 
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I. INTRODUCTION 

Periodic crystalline surfaces are a convenient testing ground for modeling surface evolu- 
tion and considerable efforts have been devoted to their study. The decay of one and two 
dimensional surface modulations has been investigated extensively. The emerging exper- 
imental picture is that above the roughening transition temperature, Tn, such structures 
decay as pure sine wavesa. This fact is well explained by Mullins' theoryH which treats the 
surface as an isotropic object. However, below Tr the situation is different and one must 
take into account the anisotropy of the underlying crystal lattice. This anisotropy produces 
cusp singularities inJhp crystal surface tension at high symmetry planes, which lead to the 

formation of facetsa'ETQ. 

From the theoretical point of view, the existence of cusps in the surface tension com- 
plicates the description of surface evolution in terms of continuous models. Such models 
usually assume that mass transfer follows gradients in the continuous surface chemical po- 
tential. However, the chemical potential is singular at facet edges and the dynamics in the 
vicinity of these regions must be treated with special care. 

There are basically two ways to deal with the above problem. The first one is to avoid 
it all together_by resorting to microscopic models. Step flow modelsLTEl and Monte Carlo 
simulationalMcJ describe surface dynamics in terms of discrete microscopic objects. On this 
level the surface chemical potential is no longer continuous and there are no singularity prob- 
lems. Facets are natural phenomena in such a framework. However in going to microscopic 
models one usually loses the mathematical simplicity of the continuum approach. Step flow 
models and Monte Carlo simulations are usually solved numerically and it is difficult to 
derive the large scale dynamics from them. 

The_second way is ±o use continuum models in which singular points are treated as 
movingEj or stationaryQ boundaries. Alternatively, surface tension singularities can be 
smoothed outE-O'ta assuming that if smoothing is done on a small length scale the large 
scale results will not be affected. The drawback of these continuum approaches is that by 
starting with a continuum rnjpdel one loses sight of the underlying microscopic kinetics. As 
pointed out by Chame et al.a, the connection of these models with the details of the micro- 
scopic processes is unclear. We show below that these details are important and may lead 
to large scale effects. 

In this work we re-examine the relaxation process of unidirectional surface modulations. 
Our aim is to derive a continuum model for the surface evolution which is consistent with 
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the microscopic kinetics. We adopt here the same approach we used in the study of other 
surface structures^rE3. First we study the system's behavior in terms of a standard step 
flow model. The model allows steps of opposite sign to interact; such interactions can 
lead to facet formation. Numerical simulations of the step model suggest that the surface 
slope D(x,t) follows the simple scaling law D(x,t) = a(t)F(x). The derivation of the step 
model and its numerical analysis are carried out in Section II. In Section III we use the 
scaling ansatz to transform the step flow model into a continuum model, i.e., we derive 
the differential equations for the functions a(t) and F(x) directly from the step velocities 
of the discrete model. In Section IV we solve the resulting differential equations and find 
the scaling function F. We find an impressive agreement between this scaling function and 
simulations of the step model. 

In Section V we discuss the relevance of our model to AFM measurements of decaying 
surface modulations taken by Tanaka et alM We show that the best fit solution of our model 
agrees with the experimental data. On this basis we derive a formula which connects micro- 
scopic parameters in the system with the measured life time of the profile. Our conclusions 
and their relation to existing work are presented in section VI. 

II. STEP FLOW MODEL OF ID MODULATIONS 

Consider a ID grating of wavelength A fabricated on a high symmetry crystal orientation 
as illustrated in Fig. |l|. Below the roughening transition one can ignore the existence of 
islands and voids and describe the surface as consisting of flat terraces separated by straight 
atomic steps. In the absence of evaporation and when no new material is added, these 
steps move by mass exchange with the adatom diffusion fields on neighboring terraces. Such 
evolution was treated in the spirit of the Burton-Cabrera-Frank modeled by a number of 
authors. In what follows we use a standard step model for surface evolution. This model 
is essentially equivalent to the one used by Ozdemir and Zangwill in Ref. [7] except that we 
allow steps of opposite sign to interact. For completeness we now derive the model. 
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FIG. 1. Illustration of a modulated surface below the roughening transition. The step labeling 
is indicated. We label the terraces such that the nth terrace separates the nth and n + 1th steps. 

To describe step motion mathematically, one has to solve the adatom diffusion problem 
on the terraces. Let us denote the adatom concentration field on the nth terrace by C n {x). 
In most situations, the time scale associated with step motion is much larger than the time 
scale of surface diffusion. One can therefore assume that the adatom concentration field is 
always in a steady state; i.e, for any given step configuration, C n (x) reaches a steady state 
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before the steps move significantly. Within this quasi-static approximation C n (x) obeys the 
static diffusion equation V 2 C n (x) = 0, which is readily solved by 

C n (x) = a n + b n x . (I) 

Next we look at the boundary conditions for these diffusion fields. Near the step edges 
the diffusion currents are determined by the flux of atoms emitted or absorbed by the steps. 
Assuming linear kinetics characterized by an a^chment- detachment rate k, the current of 
atoms on both sides of the nth step is given bj 



J2-U 



J n _i = -L> s VC n _ a (x n ) = k [C„_! (i n ) - C e n q ] , 
J n = -D s VC n (x n ) = -k [C n (x n ) - C%\ . (2) 

Here D s is the adatom diffusion coefficient, x n is the position of the nth step, J n is the 
current across the nth terrace and C^ q is the equilibrium adatom concentration at x n . 

Eqs. (HI) and @ can be used to calculate the constants a n and b n and thus fully determine 
the adatom concentration fields. We find that 

Jn = D s b n = ~2£)~ s ' • (3) 

i. i -En+l Xn 

Mass conservation implies that the velocity of the nth step takes the form 

<l> '" Q(Jn-Jn-l), (4) 



dt 

where Q is the atomic area of the solid. 

In order to complete the model, we write down expressions for C^ q . To this end we 

introduce the step chemical potential fi n , which is associated with the addition of an atom 
to the solid at the nth step. C^ g depends on the step chemical potential in the following 

way: 

C ^e*p^=^(l + ^). (5) 

C eq is the equilibrium concentration of a noninteracting step, ks is the Boltzmann constant 
and T is the temperature. 

Finally, to evaluate fi n we take into account repulsive interactions with strength (3/2 
between nearest neighbor steps of the same sign: 

U(x n ,x n+1 ) = — -2- (6) 

^ (^n+l — x n ) 

Such a repulsion is consistent with entropic as well as elastidl@ interactions between straight 
steps. We also take into account the possibility of an attractive interaction with strength 
(3/2 between nearest neighbor steps of opposite signs: 

U {x n , x n+1 ) = -— ^— . (7) 



An elastic attraction of this form may exist in some ma±erialsc3. In addition, step-antistep 
annihilation events are accelerated by step fluctuationsBtE In our one dimensional model, 
step fluctuations are not treated. Instead, this kinetic effect can be introduced as an effective 
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step-antistep attraction. Although the form of this kinetic attraction may be different, Eq. 
(0) can be viewed as a phenomenological ingredient of our model which enables us to study 
the effect of step-antistep attraction on surface evolution. Using the above interactions the 
chemical potentials of steps 1, . . . , N in Fig. [I] are given by 



d\U(xo,x 1 ) + U(x 1 ,x 2 )\ (3 (3 



+ 



dx 1 (x 2 - xi) 3 (x 1 - x ) 3 

d[U (x n _ u x n ) + U ( 

)\ P $ o v , 

9x n {Xn+l ~ x n) (%n ~ ^n-l) 



d 

flN = — 



U(x N -x,x N ) + U(x N ,x N+ i) p p 



dx N ( Xn - xat_i) (x n+1 - x N 



,3 - 



(8) 



Eqs. @, (f|), (H) and (|) constitute a complete model for surface evolution. This model 
depends on six parameters but many of them are trivial. The only non trivial parameters in 
the model are the ratio g = (3/(3, which measures the relative strength of the two interactions, 
and the length I = The latter determines the rate limiting process in the system. When 
/ — > attachment-detachment events are relatively fast and the kinetics is diffusion limited 
(DL). In the opposite case, / — > oo, diffusion is relatively fast and the kinetics is attachment- 
detachment limited (ADL). The combined effect of all the other parameters can be scaled 
out by transforming to dimensionless variables. In the DL case we choose 

x - ^x 
X 

~ ( X \ 5 k B T 

t = — ^ 1 . 9 

\2nJ QD s C e ip 

In all other non-diffusion-limited (NDL) cases we choose 

t (-) -J^—t. (10) 

with the same definition of a;. In the remainder of this section we study surface evolution 
through numerical simulations of the above model in terms of the dimensionless variables x 
and i. 

Figure g presents the simulation results of a typical system. We show a plot of the step 
configuration in one wave length of the profile. The initial step configuration corresponds 
to a sinusoidal surface profile. After a short transient step-antistep pairs at the top and 
bottom terraces start annihilating. As a result, the steps in the sloping parts of the profile 
become less densed. 
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FIG. 2. Simulation result of a typical step system. Solid (dashed) curves show the position of 
steps (antisteps). Steps and antisteps collide and annihilate at the profile maximum (x = vr/2) and 
minimum (x = —ir/2). 

To further understand the surface evolution let us study step kinetics through their 
density function (the profile slope). We define the step density between two neighboring 
steps of the same sign as the inverse step separation, i.e., 



D f x n+1 +x n - 



(t) - x n (t) ' 



(11) 



where the ± signs distinguish between steps and antisteps. The density between two steps 
of opposite sign is defined to be zero, consistently with the profile slope. 

D(x,t) is a continuous function of t, defined at a discrete set of positions at any given 
time. It turns out that the evolution of this function obeys a simple scaling ansatz. When 
we scale the step density to unity at its peak and then plot density functions from different 
times, all the data collapses onto a single curve. In other words, there exist functions F(x) 
and a(i) = max^} D(x, t), which satisfy the relation 



D (5, t) = a (t) F 



x 



(12) 



for every x = (x n+ \ + x n ) /2. Note that since the step positions vary continuously with time, 

F is a function of a continuous variable. 

Figure f| demonstrates the data collapse obtained when the density function is divided 

by its amplitude, ait), at different times for the ADL case with (Fig. [5]a) and without (Fig. 
|||b) step-antistep attraction. Only half a period of a profile modulation is shown, since 
by symmetry the density of antisteps in the other half behaves in exactly the same way. A 
similar behavior is seen in the DL limit. We note that in all cases the system exhibits scaling 
to a good approximation, but the quality of the data collapse is slightly better in the cases 
without step-antistep attraction (g — 0). The second observation is that when step-antistep 
attraction is present, the scaling function is steeper and has a smaller step density near the 
profile extrema. This is a result of the faster step-antistep annihilation process due to the 
attraction. 
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FIG. 3. Data collapse of density functions in the ADL case with a) g = and b) g = 24. Note 
that in the presence of step-antistep attraction (g = 24) the density function is steeper and has a 
smaller value near the profile extrema. 

While there is no significant difference between the ADL and DL scaling functions, the 
time dependent amplitude, az(t), does show different behaviors. In Fig. |] we plot a(t) for 
ADL and DL systems, with and without step-antistep attraction. In the ADL cases, a(t) 
can be fitted by an exponential decay, while in the DL cases a(t) ~ (t — to)" 1 consistently 
with Ref. 0- 



a) b) 
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FIG. 4. The amplitudes a(t) in different systems (circles): a) ADL with g = and g = 24, b) 
DL with g = and g = 24. In the ADL cases the amplitudes can be fitted (solid lines in a) by 

a = A* e~ l l T with r = 1.42 for the g = case and r = 1.196 for the g = 24 case. In the DL cases 
the amplitudes can be fitted (solid lines in b) by a = r/(t — to) with r = 0.785 for the g = case 
and r = 0.643 for the g = 24 case. Note that in the presence of step-antistep attraction the decay 
is faster. 

III. SCALING ANALYSIS AND A CONTINUUM MODEL 

The simulation results from the previous section suggest that the step density in the 
discrete step model can be described by the continuous functions a(t) and F(x). This 
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situation is similar in many aspects to th^oije encountered in a different step system which 



describes the decay of a crystalline coneE3E!l. Although the cone system obeys a different 



scaling scenario, we use here the same scaling analysis, described in detail in Ref. [21]. Our 
aim is to derive the differential equations for the functions a(t) and F(x). 

Assuming that the scaling ansatz ( |I2"D indeed holds, we write the full time derivative of 
the step density: 

dD _dD 3D dx 

~dT ~ ~dt + ~b~x~ ' H ' ( ' 

Equation (|TJ|) can be evaluated in the middle of the terrace between two steps, i.e. at 
x = (x n + x n+ i)/2. The time derivative of x is then given by dx/dt = (x n + x n+ i)/2 where 
x n = dx n jdt. The l.h.s. of Eq. flT3| ) can be calculated by taking the time derivative of Eq. 
(|TT|): dD/dt = —D 2 {x n +i — x^j. Using the scaling ansatz, flT2|), we express Eq. flT3| ) in 
terms of the functions a(t) and F(x): 

aF , x n+1 +x n + ^ F + a2p2 ^ _ ~ n ) = o . (14) 

a and F' are the t and x derivatives of a and F, respectively. The step velocities x n and 
x n+ i can in principle be expressed in terms of the x„'s using Eq. @, but we defer this 
manipulation to a later stage. 

Let us also rewrite Eq. ([□]) in terms of a and F: 

X n+l ~ x n = r r/ ~ : ~ \ TTvi • l-'-^'J 

aF[(x n+ i + x n )/2J 

According to this, the difference between successive x n 's is of order a~ l wherever F does 
not vanish. This allows us to expand Eq. ( |14"D in powers of a -1 when a is large and the step 
density is high. We can consider only the leading order in a -1 , an approximation which 
becomes exact in the scaling limit a — > oo. The differences x n +k ~ x, where x denotes 
the middle of the terrace, are also small as long as k is finite. This allows us to go to a 
continuum limit in the variable x in the following way. We evaluate the function F at the 
position (x n+ k + x n+ k+i)/2 by using its Taylor expansion 

Xn+k + Xn+k+l \ a 1 1 d m F(x) ( X n +k + X n +k+l -A™ / 1 „\ 

~ ~ = —\ A~m O X l ■ ( 16 ) 

Xn+k+1 ~ Xn+k m= Q TTV. dX \ A / 



Now we expand 



oo 

X n+k = X + <t>km®~~ m , (17) 
m=l 



and insert the expansion into Eq. (|TJJ). By requiring all orders of a' 1 to cancel, the coeffi- 
cients (pkm are found for any values of k and m. These coefficients involve the function F 
and its derivatives evaluated at x — (x n + x n+ i)/2. 

We now return to Eq. ([b|). It depends on the velocities x n and x n+ ±, which in turn 
depend on x n -2, ■ ■ ■ ,x n +3- Using Eq. (|T7|) we expand Eq. (0) in a^ 1 . The result of this 
expansion in the general NDL case is 
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' 1 X -".F-O(a-*)-0. (18) 



dx 2 J a 

In the DL case the expansion result is different: 

3 d 4 F 2 a „ m / ,\ . i 

2'lar + ^' F + ( a )= ' < 19 » 

Consider Eq. in the scaling limit, a — > oo. The term can be ignored and we 
are left with the first two terms, which must cancel each other. This implies that a/a is 
independent of time, i.e., in the NDL case 

a oc e ~ i/TNDL , (20) 

where tjvdl is the life time of the profile. Note that the scaling limit a — > oo corresponds 
to very early times t — > — oo. The equation for the NDL scaling function in this limit is 




d 2 F 2 



dx 2 I T NDL 



. (21) 



At any finite time when a itself is finite, there will be a correction to scaling of order a~ 2 . 

The same arguments can be applied to Eq. ([Tj]) in the DL case. We find that a /a 2 is 
independent of time, i.e., 



a = . (22) 
t-t 

In this case the scaling limit corresponds to a finite time, t — > tj. The differential equation 
for the DL case is 

3 d 4 F 2 F , . 

r-« = • 23 

2 dx* r DL 1 ; 

Again there is a correction to scaling of order a~ 2 when t > to- 

Note that the scaling equation is the same for all NDL cases. Thus, for all non-zero 
values of I = the scaling behavior is identical to that of the ADL limit. Since in the 
NDL case the terrace sizes vanish in the scaling limit, diffusion across terraces is fast and, 
unless the system if purely DL (/ = 0), attachment-detachment is always the rate limiting 
process. This fact makes the DL case special and somewhat delicate. Since in real systems 
the diffusion coefficient is always finite, the scaling limit of real systems will show ADL 
behavior. Nevertheless, at finite times we can expect crossover from ADL to DL behavior. 
This should happen at the time when « I. At later times, when 3> /, the system 
should show DL behavior and obey Eq. (|23"D , provided a is still large enough to neglect the 
a~ 2 corrections to scaling. 

At this point we have a continuum model derived directly from the discrete step equa- 
tions of motion. In the scaling limit our model provides an exact connection between the 
microscopic step kinetics and the macroscopic surface evolution. It is interesting to compare 
it with other continuum models, which do not emerge from the underlying step model. 

In the usual continuum approach one assumes that surface dynamics is driven by the 
surface tendency to decrease its continuous surfaceJtension. Below the roughening transition 
this surface tension has a cusp and takes the formEEZI 

a(D) = a + <t 1 \D\ + a 3 \D\ 3 , (24) 
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with D denoting the surface slope. The surface chemical potential is consequently given by 



"="^ i ' (25) 

where a' is the derivative of a with respect to the slope D. This chemiml^otential gives rise 
to surface currents, which are proportional to dfi/dx in the DL caseQ'tDEH and to d/j/dh = 
jj!^ in the ADL case@. The equation of motion for the profile slope in regions where D > 
is then given by 

dD d 2 (ld*D*\ 

K ~ 3a ^ T^^OT (26) 



dt dx 2 V .D <9x 2 



in the ADL case, and 



dD d*D 2 

-m K - 3cT3 ^ (27) 

in the DL case. 

To compare our model with the above equations, let us find the general equation of 
motion for the density function; i.e., the partial differential equation for the step density 
that by separation of variables, D(x, t) = a(t)F(x), reduces to Eqs. fl20p & (pip in the NDL 
case and to Eqs. fl22f) & (|23|) in the DL case. In the NDL case we can replace tndl in Eq. 
(PI) by —a/a. The equation can then be written as 

3 <9 2 (I d 2 D 2 \ dID_ 

l"di 2 '{D'~di 2 ~) + ^t~ 0, (28) 

which is consistent with Eq. (|26|) and confirms Nozieres' suggestionil for ADL kinetics. Note 
that Eq. (|28|) applies only to ADL cases and not to the full range of NDL systems, because 
in the scaling limit all NDL systems obey Eq. (|2TD , which describes ADL kinetics. Since Eq. 
(p8|) was derived from Eq. (^) it corresponds to ADL kinetics. 

For the DL case we replace t dl in Eq. (|23|) by —a/ a 2 . We then have 

3 d A D 2 dD n 



consistently with E q. (|27|) and Refs. and Let us however remark that the validity 
of the general Eqs. (|28|) and (|29|) is questionable. In a non scaling case these equations have 
corrections which may be important. We showed that in the scaling state these corrections 
vanish and neglected them in moving from Eqs. (pl|) & (pi]) to (pT|) & (|23|). In a general 
scenario, neglecting these terms is an approximation that must be justified. 



IV. SOLUTION OF THE CONTINUUM MODELS 

In this section we find the scaling functions in the NDL and DL cases. Since the slope- 
up and slope-down sections are symmetric, we need to consider only half a period, i.e., the 
section x G [— 7r/2, 7r/2]. The scaling function is thus symmetric about the origin. 

We now discuss the possibility of macroscopic facets at the profile extrema. In terms 
of the discrete step model every terrace is a facet. However, in the scaling limit the step 
density diverges with a(t), and the terraces in all regions where F is finite become truly 
microscopic. In regions where the scaling function is identically zero the terrace size does 
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not necessarily vanish and macroscopic facets may form. To account for macroscopic facets 
we should, in principle, allow for regions of arbitrary size with F = 0. However, simulations 
of the discrete step model indicate that if such regions appear at all, they form around the 
profile extrema where the step density is zero by definition. To account for macroscopic 
facets, we therefore allow for the existence of two special points at ±Xf acet , which denote the 
positions of facet edges. In the regions (— ir/2, — Xf ace t) and (x/ acet , 7r/2) the scaling function 
F vanishes identically. In this notation we can also describe systems without facets simply 
by setting Xf acet = vr/2. Note that the scaling function should satisfy Eq. ( f2~T|) or ( |23|) only 
in the interval (—Xf aC et,x facet) and not on the facet, since these equations are not valid when 
F = 0. 

Next we show that if there is a macroscopic facet, F is continuous across the facet edge. 
To see this, consider the annihilation process of the first step. Its velocity satisfies 

~ (x 3 - x 2 )~ 3 - 2(x 2 - x x )' 3 - g(x x - x )~ 3 

xi oc ^f— ~ ■ ( 30 ) 

*f - + x 2 - xi 

The denominator in the above expression is always positive. Therefore the direction of 
motion of the first step is given by the sign of the numerator, which must be negative 
because the first step is moving towards annihilation with the zeroth step. Assume now that 
in the scaling limit we have a macroscopic facet. This means that the distance X\ — Xq is 
finite when the annihilation process starts. Assume also that Ff ace t = F(xfacet) is finite. 
This makes the distance x% — x 2 microscopically small, since in the scaling limit it vanishes 
as 1/ (aF facet)- We can therefore neglect the last term in the numerator of Eq. fl3"0|). It is 
easy to see now that if x 2 — x\ < 2 1 / 3 / (aFf aC et) the velocity of the first step is positive. Thus 
if Ffacet is finite, the first step is bounded to the second one at a distance which is much less 
than the facet size. In such a situation the first step cannot annihilate. We conclude that 
having a macroscopic facet {x facet < 7r /2) requires the scaling function to vanish at the facet 
edge. Thus the scaling function is continuous at Xf acet with F(xf acet ) = 0. Furthermore, by 
expanding the scaling function in the vicinity of Xf ace t, it can be shown that both the NDL 



and DL scaling functions can be expressed as power series of yjx facet — x: 

~ / r- 

F (x) = J2 a n X facet ~ X) . (31) 

71=1 

Hence with a no n- vanishing coefficient ai, the derivatives of the scaling function diverge at 

Xfacet- 

Next we specify boundary conditions for the scaling function F . Since Eqs. and (p3|) 
are fourth order differential equations we need four boundary conditions in order to solve 
them. Three conditions are set by the normalization and symmetry of the scaling function: 

F(0) = 1 , 
F'(0) = , 

F"'(0) = . (32) 

A fourth boundary condition can be found by considering mass transfer in the decaying 
profile. Mass is transfered from the decaying peaks to the valleys through the sloping parts 
of the profile. The step-antistep symmetry of the system excludes mass transfer through 
the profile extrema. We can thus calculate the (dimensionless) flux through the origin by 
calculating the (dimensionless) volume change in the interval [0,7r/2]: 

d W 2 

J(Q)= h{x,t)dx . (33) 



dt Jo 
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h(x, t) is the surface profile measured in steps. In the continuum limit it is given by 



h(x, t) 



(34) 



Integrating Eq. ( j33j) by parts and using the scaling form (|T2| ) together with the possible 
existence of zero density facets results in 



J(0) 

To evaluate the integral in Eq. 
and integrate from zero to Xf acet : 



a 



■■facet 



F ■ (tt/2 - x ) dx . (35) 
in the NDL case, we multiply Eq. (PH) by (tt/2 — x) 



^ F • (tt/2 -x)dx = ^ 



(36) 



The r.h.s. of the last equation can be integrated by parts. Inserting the result into Eq. (|35] 
we find that 



J(0) 



3r NDL a 



(tt/2 - x) 



d_ 

dx 



'1_ d 2 F 2 " 
F 



dx 



+ 



1 d 2 F 2 



F di 2 



j facet 



(37) 



Another method for calculating the flux through the origin relies on the discrete step 
system. The adatom flux across the nth terrace is given by Eq. (|3|), which can be expanded 
in powers of a~ l using Eq. (|T?|) . In our dimensionless variables, the leading order of this 
expansion in the NDL case is 



J(x) 



3a d 2 F 2 
AF 



dx 



(38) 



Evaluating this expression at x = and comparing with Eq. fl3"7p, we obtain the boundary 
condition 



(tt/2 -x) 



d_ 

dx 



\_ d 2 F 2 " 
F 



dx 



+ 



1 d 2 F 2 



F di 2 



. 



(39) 



facet 



In deriving this boundary condition we used the fact that F'(0) = F"'(Q) = and that in 
the NDL case t^dl = — a/a. 

The DL case is very similar. We can evaluate the flux through the origin both from Eq. 
(j35|) and from a direct expansion of Eq. (^|). Equating the two results we obtain the fourth 
boundary condition in the DL case 



(tt/2-x) 



d 3 F 2 d 2 F 2 
+ 



dx 3 dx 2 







(40) 



facet 



At this point we have four boundary conditions for the scaling function, three at the 
origin and one at Xf acet . In addition we know that F(xf acet ) = if Xf acet < tt/2. We may 
now find unique scaling solutions if we know the values of t^dl and t^l, which enter in the 



differential equations ( PTj ) and 
show that these life times are re 



_ . What determines the values of tnjjl and tjjl? We now 
ated to the behavior of the top and bottom steps, which are 



unique in the sense that they each have a neighboring step of opposite sign. Our continuum 
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model, which treats all steps on equal footing, does not contain any information about this 
unique behavior and therefore t ndl and t dl must be calculated or measured directly from 
the discrete step system. 

To relate the profile life times to the behavior of the top and bottom steps, consider the 
number of steps in a half slope up region of the profile. One can show from Eq. (|34|) that 
it is given by 



N(t) 



tt/2 



D(x, t)dx 



f-E facet 

a / Fdx 
Jo 



The integral I = £ facet Fdx can be evaluated by integrating Eqs. (^1|) or 
or the DL cases, respectively. The results are 



(41) 

in the NDL 



NDL 



1, 



3t ndl d ( 1 d 2 F 2 " 
4 dx \ F dx, 

Sr DL d 3 F 2 



-•facet 



DL 



dx 3 



(42) 



facet 



Combining Eqs. ([41]) and ( |4*2| ) we can calculate At, the time interval in which the system 
looses one step, i.e., the annihilation time of the top step: 



dt 1 
At = — = — . 

dN al 



(43) 



Eq. (^) relates the profile decay rate, a, to the annihilation time of a single step at the 
profile peak, At. We could have used this relation to set an additional condition at Xf acet if 
we knew At. Such a condition would select a single solution with a specific life time from 
the one dimensional family of possible solutions. However, as we mentioned above, At must 
be obtained directly from the discrete system, because our continuum model ignores the 
unique behavior of the top step. 

Unable to calculate the life times of the profiles analytically, we measured t^dl and tdl 
from our simulations. Using these values we solved Eqs. ([TO) and ( p3|) numerically. In Fig. |5| 
we show the results for a single slope-up region of the profile in the ADL case and compare 
them with the simulation data. The slope-down regions of the profile behave in exactly the 
same way, but with a negative step density (to be consistent with the profile slope). In the 
insets we show the surface profile, which is obtained by integrating the scaling function over 
alternating regions of steps and antisteps. Similar results were obtained in the DL case. 

The following picture emerges. Both in the ADL and DL cases, when step-antistep 
attraction is absent (g = 0) the unique scaling solutions, which satisfy all the boundary 
conditions, have no facets. The scaling functions obey Eqs. (|2T| ) and ( p3[ ) in the entire 
interval (— tt/2, tt/2) and are finite at the profile extrema (x = ±tt/2). This results in cusp- 
like peaks and valleys in the surface profile. The agreement between the scaling functions 
and the simulation data is excellent in these cases. When step-antistep attraction is present 
(g = 24) the situation is different. The unique scaling solutions, which satisfy all the 
boundary conditions, have small facets near the profile extrema. In the facet region the 
scaling functions are identically zero and as a result the profiles have flat peaks and valleys. 
Here we also have good agreement with simulation data, but not as good as in the g = 
cases. This is because the scaling ansatz is a better approximation when g = 0. 
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FIG. 5. Scaling function solutions (solid) compared with simulation data (circles) in the ADL 
case. Only a single slope up region is shown. The insets show the resulting surface profiles obtained 
by integrating the scaling functions over alternating regions of steps and antisteps with a) g = 
and b) g = 24. 



V. CONNECTION WITH EXPERIMENT 

In this section we test our model against experimental results. Specifically we compare 
our predictions with AFM measurements carried out by Tanaka et al@. They measured 
the amplitude of decaying ID Silicon gratings fabricated on a surface vicinal to the Si(OOl) 
plane. They report an exponential decay of the surface amplitude with life times scaling 
as the fourth power of the grating wavelength. The observed profiles had flat peaks and 
valleys. These results qualitatively resemble the behavior of our model in the NDL case in 
the presence of step-antistep attraction. 

To prepare the ground for a more quantitative comparison, let us first discuss the rel- 
evance of our model to the experimental system. The experimental conditions where such 
that the Si surface was below its roughening transition. It is therefore reasonable to expect 
step-dominated surface kinetics. However, the step structure in the experiment is different 
from the idealized one dimensional array of straight steps we considered in our model. Since 
the grating was fabricated on a vicinal surface, with the grooves direction perpendicular to 
the unperturbed steps, the steps are not straight. This can be easily shown by writing the 
surface profile as 

Z(x, y) = h(x) + ay . (44) 

Here h(x) is the one dimensional periodic profile along the x axis and Z(x, y) is the two 
dimensional surface profile. The parameter a measures the slope of the original vicinal 
surface with respect to the high symmetry (x, y) plane. The constant height contours which 
define the steps are given by 

V«{x) = ■ (45) 

a 

Evidently, we deal here with curved steps, aJ^ct which complicates surface dynamics. This 
situation was studied by Bonzel and Mullinsc3 who investigated the smoothing of perturbed 
vicinal surfaces. In addition to step-step interactions, step motion is also driven by the steps 
line tension, Y. This is reflected by the addition of a term fi™ rv (x) = T/R n (x) to the step 
chemical potential, where R n {x) is the local radius of curvature of the nth step. 
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The above arguments suggest that our straight steps model cannot describe surface 
dynamics in regions where line tension is important. To justify the application of our model 
to the experimental system we must therefore show that we can neglect the effect of step 
line tension. We do this by considering the scaling limit of the system. We show that in this 
limit the steps in the sloping parts of the profile are effectively straight. The only regions 
where step curvature is significant are the profile extrema, where our continuum model is 
not valid anyway. Unlike Bonzel and Mullins we consider situations where step curvature is 
a driving force for step straightening only at the profile extrema. 

Assume now that the step system ( fI5D exhibits scaling and obeys Eq. (0). D(x,t) = 
dh/dx is the step density along the x axis. Under these assumptions we can express the step 
curvature contribution to the chemical potential as 

Fy'n FaF' 

= r " ]3 /2 = J a2p2 ^/2 > ( 46 ) 



with primes denoting derivatives with respect to x. Note that according to Eq. (^6|), 
fj,^ trv (x) = fi curv (x) is independent of n. 

Let us estimate the effect of fi curv (x) on step kinetics. This chemical potential gives rise 
to surface currents parallel and perpendicular to the steps. The divergence of these currents 
contributes to the step velocities. Along a step we can estimate the resulting current as 
J|i oc dn curv (x) / ds where s is the step arclength. The contribution of this current to the step 
velocity is proportional to 

d 2 fi curv (x) 
cb 2 ~ 

/ aF m a 3 (3 F' 3 +10 FF' F" -2 F 2 F'") a 5 (-15 F 2 F' 3 +10F 3 F' F"-F 4 F" 



(47) 



In the scaling limit a — > oo, this contribution decays as a 4 wherever F is finite. This result 
should be compared with the step velocities of the original model. Using the expansion 
Eq. (|TT|), we can expand the original step velocities, Eq. (|), in a" 1 . The outcome of this 
manipulation is a velocity of order 1. Thus, in the scaling limit, the curvature driven current 
along the steps results in a negligible contribution to the step velocities. 

The current perpendicular to a step is proportional to the chemical potential difference 
between two adjacent steps. This can be estimated from J_l oc h ■ Vfi curv (x)/\ VZ\ where 
n is a unit vector perpendicular to the step. However, since we are only interested here 
in the leading order in a -1 , we can make the approximation h ~ x, \VZ\ rs aF and 

J± ~y dti d wherever F does not vanish. The contribution of these currents to the step 
velocity is consequently proportional to the current difference between two adjacent terraces 
which we estimate as 



1 d 1 d^ curv {x)\ ra" 4 a 2 (15F' 3 — lOFF'F" + F 2 F 



aF dx \ aF dx I F 7 



+ 0(a' 6 ) . (48) 



Again this contribution in negligible in the scaling limit. 

The above argument relies on the fact that in the sloping parts of the profile F is finite. 

This is of course not true near the profile extrema. Step curvature is therefore important 

at the facets near the profile peaks and valleys where the scaling function F is identically 

zero. Thus, we can divide the system into two types of regions: the sloping parts where our 

continuum model is valid and the facets where our continuum model breaks down. Recall 
however that our model breaks down on facets in any case, and we solve the model only on 
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the sloping parts of the profile. The kinetics of the annihilating steps on the facets enters 
the continuum model only through the life time of the profile, which we have to take from 
the discrete step system. We can use this strategy here, i.e., by measuring the life time 
of the profile from the experimental results we can predict (using Eq. the density 

scaling function. Alternatively, by fitting the scaling function to the experiment we can 
make a connection between the experimental life time T exp , and the microscopic parameters 
of the system. In Fig. ^| we show a comparison between a snapshot of the normalized 
experimental slope and the best fit solution of Eq. (|2l|) . This solution corresponds to a 
life time t ndl = 0.945 in dimensionless units. According to Eq. fllCf) the three microscopic 
parameters k, C eq and (3 satisfy 

4 k B Tr NDL ^ (4g) 

The parameter values for the experimental data of Fig. ^| are A = 5 fim, r exp ~ 2.4 • 10 3 
minutes and T = 900°C. 





-2tc 2n 

X 

FIG. 6. Comparison between a snapshot of the experimental slope and the best fit solution of 
Eq. (]2l|). This solution corresponds to a life time tndl = 0.945. The experimental slope was 
normalized to a 2ir wavelength and unit maximal slope. 

Although Fig. |H shows good agreement between theory and experiment it does not by 
itself provide enough evidence for scaling in the experimental system. To verify that the 
experimental system scales in time one must analyze several snapshots of the surface slope 
taken at different times. However, the fit quality and the fact that the surface amplitude 
decays exponentially with a life time proportional to the fourth power of the grating wave- 
length, strongly suggest scaling. 



VI. SUMMARY AND DISCUSSION 

In this work we have studied the relaxation process of one dimensional surface mod- 
ulations below the roughening transition temperature. Simulations of a step flow model 
suggest that the discrete step density function, D(x,t), exhibits scaling throughout most of 
the decay process, for a wide range of the model parameters, spanning DL to ADL kinetics 
and a range of strengths of the attraction between steps of opposite sign. 

Relying on the above observations, we then transformed the discrete step model into a 
continuum model for surface dynamics. This was done using the scaling ansatz D(x,t) = 
a(t)F(x), and we were able to write down the differential equations for the functions a(t) 
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and F(x). These equations were derived directly from the discrete step equations of motion. 
We showed that they are exact in the scaling limit. 

We investigated the boundary conditions for the scaling function F and found F numer- 
ically. A comparison between the resulting solutions and density functions from simulations 
of the discrete step model shows impressive agreement. 

u Finally, we applied our continuum model to experimental data measured by Tanaka et 
al@. They reported exponential amplitude decay of Si(OOl) gratings and profiles with flat 
peaks and valleys. In their experiment the profile life times scaled as the fourth power of 
the grating wavelength. These properties resemble the behavior of our model in the NDL 
case with step-antistep attraction. We showed that in the scaling limit, our one dimensional 
model is adequate in spite of the fact that the steps in the experimental system are not 
straight. We compared the best fit solution of our model to a normalized snapshot of the 
experimental profile slope and found good agreement. On the basis of this agreement we 
suggested that the profile exhibits scaling and made a connection between the measured 
profile life time and microscopic system parameters. 

We now discuss the main aspects of our results and put them in perspective with existing 
work. In the scaling limit our continuum model is an exact result of step kinetics. We showed 
that the resulting general differential equations and (P§D) are equivalent to existing 

phenomenological modelsB'utH, thereby confirming that these models are consistent with 
step flow. Although we have not discussed this issue here, it can be shown that the scaling 
solutions D(x, t) = a(t)F(x) are linearly stable solutions of the general differential equations 
((|28|) and ([29])). This fact may explain why scaling solutions are selected. However, our 
analysis shows that these differential equations have corrections which originate from the 
discrete nature of the system. When these corrections are important, one must take into 
account the full expansion series (0). This results in differential equations of infinite order. 
Therefore Eqs. (|28|) & (p9|) and other equivalent models are valid only in situations where 
corrections can be ignored. We proved that the scaling scenario is such a case by showing 
how these corrections vanish in the scaling limit. 

The second point we emphasize is that even in a scaling limit there are regions where 
the discrete nature of steps becomes important. Such regions are the profile extrema where 
step-antistep pairs annihilate. The reason for this is twofold. First, the steps near the profile 
extrema are unique in the sense that they have a neighboring step of opposite sign. They 
therefore follow unique equations of motion. This does not enter in our continuum model 
which treats all steps on equal footing. Secondly, our continuum model breaks down near 
zeros of the scaling function. This can be seen in the derivation of the model, which relies 
on F being finite. The kinetics of steps in regions where the step density vanishes are not 
described by the model, since corrections to the scaling function become important there. 
In cases where the profile extrema are faceted, our continuum model cannot account for the 
kinetics of steps on the facets even if these steps follow the general discrete equations of 
motion. 

It is therefore obvious that our continuum model (and all other equivalent differential 
equations) does not constitute a complete description of surface evolution. A manifestation 
of this statement is the fact that the relative attraction strength, g, is not a parameter of the 
continuum model, even though it strongly affects surface evolution. Since our continuum 
model breaks down near the profile extrema, their effect must be taken into account by 
supplying additional information. In our case, this information is the value of the profile life 
time. As we showed in section IV, this life time determines the annihilation rate of step- 
antistep pairs and can be considered as the last boundary condition required for finding 
the scaling function. Thus, in a scaling scenario, the effect of step kinetics near the profile 
extrema reduces to a boundary condition for the scaling function at ±Xf acet . In order to 
solve for F consistently with the discrete model, it is sufficient to supply the profile life time. 
We did this by measuring the life times of the discrete systems. 

In this respect our work is different from other continuum treatments of modulated sur- 
faces. As pointed out by Chame et aE3, lack of consistency with the kinetics of facet steps is 
the weakness of existing continuum models. One remarkable outcome of these inconsisten- 
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cies is thg controversy regarding. the appearance of facets at the profile extrema. Ozdemir & 
ZangwillQ and Hager & SflohnllS both used continuum models, which are equivalent to Eq. 

Bonzel and Preussta use a slightly different model in which the surface currents are 
proportional to the derivative of the chemical potential with respect to the surface arclength. 
But they all start from a continuum approach and do not specify the step kinetics on the 
facets. 

Ozdemir and Zangwill assume that Eq. (|29"D is valid in the full interval between the profile 
extrema and thus exclude the possibility of facet formation. Their results show smooth (but 
nonanalytic) profiles. This corresponds to some weak attractive step-antistep interaction 
(in the absence of such interaction we get cusped profiles). They showed that Eq. ( |29| ) 
admits shape preserving solutions, equivalent to our scaling solutions in the DL case. In 
particular they showed that the amplitude a of a shape preserving profile obeys Eq. ([22|). 
Our results for the decay rate in the DL case are mathematically identical, but we interpret 
them differently. While Ozdemir and Zangwill look at the long time limit in which a ~ £ _1 , 
we recognize that Eq. (|29|) is valid only in the scaling limit, where t •> t . In this case 

a ^ r 1 . 

Hager and Spohn, on the other hand, allow for facet formation by solving Eq. (|29|) with 
moving boundaries. They assume that the chemical potential, the current and the current 
divergence are all continuous at the facet edge. In their model, the chemical potential at 
the facet edge is not a step chemical potential. It is a layer chemical potential which reflects 
the free energy change caused by the annihilation of the top step-antistep pair. In this way 
they allow for the top terraces to peel rapidly. Their model produces facets which grow in 
time. 

Our results are different. First, since we start from step flow, the chemical potential in 
our continuum model is a step chemical potential. On a facet our chemical potential is not 
defined (and therefore not continuous) since there are no steps there. One can argue that 
the adatom chemical potential on the facet is equal to the chemical potential of the top first 
step. But this first step is unique and is not treated by the continuum model. It is also 
far from equilibrium with its neighbors so there is no reason to assume continuity of the 
chemical potential. Second, the appearance of facets in our model depends on strength of 
the step-antistep attraction. 

Bonzel and Preuss deal with facets differently. They round the cusp singularity of the 
surface tension, Eq. (0), and apply the kinetic surface equation everywhere. They thus 
obtain analytic profiles with very flat extrema but without actual faceting. Again the weak- 
ness of this approach is that it is not clear how this rounding is related to the microscopic 
behavior of steps. 

We conclude that our continuum model is fully consistent with step kinetics, but is 
restricted to scaling scenarios. It may serve both as a starting point and a test case, for 
new continuum models which will attempt to describe surface evolution in general and facet 
kinetics in particular. 
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